JOURNAL  or  GEOPHYSICAL  RESEARCH.  VOL.  95.  NO.  1)5,  PAGES  6779-6786,  MAY  10,  1990 


AD-A228  685 


Obtaining  Earth  Surface  and  Spatial  Deflections  of  the  Vertical 
From  Free-Air  Gravity  Anomaly  and  Elevation  Data 

Without  Density  Assumptions  DTIC 

David  M.  Gleason  \  ELECTE 

Geophysics  Laboratory,  Hatiscom  Air  Force  Base,  Bedford,  Massachusetts 


Moritz  (1980)  presents  a  density-free  scheme  allowing  for  the  analytical  or  regular  continuation  of  a 
given  set  of  free-air  gravity  anomalies,  referenced  to  the  ground,  to  any  desired  level  surface  if  a 
corresponding  set  of  elevations  (e.g.,  above  mean  sea  level)  is  available.  An  efficient  spectral 
implementation  of  this  scheme  is  discussed  by  Sideris  (1987).  A  subsequent  spectral  execution  of  the 
planar  Vening-Meinez  equation  on  the  continued  anomalies  yields  deflections  of  the  vertical  on  the 
chosen  level  surface.  The  deflections  are  brought  back  to  the  Earth’s  surface  via  a  spectrally 
implemented  Taylor  scries.  Deflections  at  a  constant  altitude  above  the  level  surface  are  obtained 
through  a  routine  spectral  execution  of  the  planar  upward  continuation  integral.  Two  sites,  having 
diverse  topographies,  were  surveyed  for  1  arc  min  by  1  arc  min  mean  free-air  anomaly  and  elevation 
values  and  for  smaller  sets  of  astronomically  determined  deflections  to  serve  as  control  or  “truth” 
values.  In  a  topographically  tranquil  but  gravimetrically  turbulent  Oklahoma  site  the  overall  RMS  of 
the  differences  between  true  and  predicted  deflections  was  0.3  arc  secs  and  in  a  rugged  New  Mexico 
site  it  was  0.6  arc  sec.  Accurate  first  derivative  terms  (in  both  continuation  steps)  require  a  I  arc  min 
data  set  as  interpolation-free  as  possible.  A  I  arc  min  data  grid  is  shown  to  be  insufficient  for 
meaningful  computations  of  the  higher  order  series  terms.  Potential  pitfalls  of  the  two-dimensional  fast 
Fourier  transform  pair  are  discussed  with  an  emphasis  on  unwanted  circular  convolution  effects 
which,  if  unaccounted  for,  can  increase  the  error  in  individual  predicted  deflections  by  as  much  as 
100%. 
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I.  Introductory  Theoretical  Discussion 

Runge’s  theorem  [cf.  Moritz,  1980]  states  one  can  always 
find  a  harmonic  function  T*,  arbitrarily  close  to  the  Earth’s 
external  (and  harmonic)  disturbing  potential  T,  that  can  be 
analytically  or  regularly  continued  (be  it  upward  or  down¬ 
ward)  from  the  Earth’s  actual  surface  to  any  desired  level 
surface.  Therefore  a  converging  Taylor  series  links  the 
observed  A gP  free-air  gravity  anomaly,  referenced  to  the 
ground  and  given  by 

«  observed  at  P  ~  Yat  the  tclloroidal  point  R 

to  the  anomaly  A #'/»•,  referenced  to  the  chosen  level  surface 
(see  Figure  1  for  subscript  referrals).  Heiskanen  and  Moritz 
[1967]  define  the  tclluroid  as  a  near-ground  surface  whose 
normal  potential  at  R  is  equal  to  the  actual  potential  at  the 
corresponding  ground  point  P  (the  points  P  and  R  reside  on 
the  same  ellipsoidal  normal)  and  show  how  to  compute  the 
normal  gravity  yK. 

Under  such  a  regular  continuation  of  T,  the  level  surface 
A,g’  set  reflects  the  Earth's  exterior  gravity  field.  Inserting 
the  A f;'  set  into  Stokes'  formula  yields  a  potential  7"  which  is 
harmonic  above  the  chosen  level  surface  (thus  unrelated  to 
the  Earth’s  nonharmonic  interior  field)  and  which  agrees 
with  the  actual  T  on  and  above  the  Earth’s  surface.  There¬ 
fore  any  masses  existing  outside  the  chosen  surface  are,  in 
effect,  shifted  to  its  interior  without  changing  T on  and  above 
the  ground.  This  regular  "continuation"  of  T  (and  thus  Aj t?) 
is,  at  least  theoretically,  preferable  to  traditional  Bouguer 
and  isostatic  "reductions"  because  the  former  can  be  ac¬ 
complished  without  any  assumptions  regarding  the  densities 


of  such  masses  and  without  approximating  the  vertical 
gradient  of  actual  gravity  by  its  normal  counterpart  dy/dz  a  4 
-0.3086  mGal/m.  As  shown  later  in  this  section,  the  vertical 
gradient  of  anomalous  gravity  is  actually  solved  for  in  the  * 
continuation  process.  ^ 

The  Taylor  series  of  interest  is  given  by 

3A  g'r  1  d2A  & 

AgP  =  A #?/>■  +  zP  — —  +  ~Zp  ----- 
dZ  2!  dZ" 


1  d3A  gp‘ 


H - Zp' 

3!  P 


hz 


.} 


(1) 


where 


zP  =  hp  -  /i*  (2) 

with  the  heights  hr  and  h'  shown  in  Figure  1. 

An  inversion  of  equation  (I)  is  desired,  i.e..  a  solution  for 
Ag '  in  terms  of  the  observed  Aj>.  Armed  with  a  set  of  Ag'  on 
the  level  surface  approximated  by  an  infinitely  extended 
plane,  one  can  compute  the  height  anomaly  f,  the  deflection 
of  the  vertical  component  along  the  meridian,  and  along 
the  prime  vertical,  tj,  at  any  point  P‘  located  on  the  level 
surface  via  the  planar  versions  of  the  Bruns  and  Vening- 
Meinez  equations  [cf,  Botnford,  1971] 
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Fig.  1.  The  applicable  surfaces. 
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where  y  is  some  mean  value  of  normal  gravity  (say,  979,800 
mGal).  Just  as  (I)  returns  one  to  \gP  from  A#'/..,  use  of 
similar  Taylor  series  allows  one  to  obtain  ( P ,  £P  and  r)P  from 
ip  -  £/>•  and  t\p 

Moritz  [1980]  shows  that  for  any  point  Q'  on  the  level 
surface 


Agy  -  2  Hq-  (5) 

n  -  0 

where  the  g£>-  (superscript  combinations  involving  n  or  A  do 
not  denote  exponentiation  while  those  involving  r  only  do) 
are  computed  recursively  via 
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where  the  Lr  values  are  also  obtained  recursively  via 
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for  any  arbitrary  k  and  the  function  L(  is  defined  below  by 
equation  (9).  The  starting  values  needed  to  implement  (6) 
and  (7)  are  given  by 


and 
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dx  dx 


(9) 


Notice  that  (9)  is  the  planar  version  of  the  well-known 
integral  expression  for  the  anomalous  radial  (vertical)  gradi¬ 
ent  < :tSg/i)r  (or  riSg/hz)  wherein  Ag  is  assumed  to  be  a 
harmonic  function  in  space.  Under  the  planar  approximation 


where  Ag  =  -dTIdz.  such  an  assumption  is  appropriate  |see 
Moritz,  1980,  p.  I76|.  The  geometric  kernel  of  (9).  1/ 
(distance) J,  dictates  0A glitz  to  be  a  high-frequency/ 
short-wavelength  phenomena.  Schwartz  (1984).  assuming  a 
widely  used  model  spectrum,  shows  that  over  609?  of  the 
value  of  ilSglitz  is  determined  by  the  behavior  of  the  gravity 
field  in  the  immediate  3  arc  min  area  surrounding  the 
computation  point.  Thus  the  g°  or  Ag  data  in  this  immediate 
area  must  be  highly  reliable  to  accurately  compute  AAg I  Hz. 

Now  equations  (3)  and  (4)  can  be  written  as 
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Thus  before  the  multiple  integrations  of  ( 10)  and  (II)  can  be 
executed,  multiple  integrations  corresponding  to  every  point 
Q'  on  the  level  surface  are  required  to  compute  the  input  g". 
These  practical  integration  concerns  can  be  overcome  by 
efficiently  executing  the  whole  scenario  outlined  in  this 
section  in  the  frequency  domain.  Sideris  [1987]  (also  see 
Sideris  and  Schwartz  [1986,  1988|  for  related  works)  uses 
analytically  defined  spectra,  as  well  as  discrete  two- 
dimensional  fast  Fourier  transforms  (FFTs),  to  compute  the 
convolutions  of  equations  (9M11).  This  “analytical"  ap¬ 
proach  will  be  summarized  in  the  next  section.  An  alterna¬ 
tive  approach,  suggested  by  Wang  [1988],  is  to  compute  all 
of  the  required  spectra,  including  the  geometrical  kernels, 
via  the  discrete  two-dimensional  FFT  pair.  Experiments  by 
this  author  have  shown  that  use  of  the  analytically  defined 
transfer  functions  yield  appreciably  more  accurate  deflection 
predictions  and  thus  will  be  incorporated  in  this  paper. 

Section  3  will  discuss  potential  pitfalls  in  any  use  of  the 
discrete  two-dimensional  FFT  pair,  highlighting  the  un¬ 
wanted  phenomena  known  as  circular  convolution.  Finally, 
section  4  will  examine  sets  of  predicted  deflections  of  the 
vertical  in  highly  diverse  topographic  settings.  Required 
input  for  the  predictions  are  a  truncated  set  of  geopotential 
coefficients  and  gridded  data  bases  of  mean  free-air  gravity 
anomalies  and  elevations. 


2.  Using  Analytic  ally  Defined  Spectra  in  the 
Computation  of  Equations  <9M  1 1 ) 

The  spectrum  of  the  rth  vertical  derivative  of  any  two- 
dimensional  planar  gravimetric  harmonic  function  h(.v,  y,  r 
=  const)  is  given  by  [cf.  Sideris,  1987] 

f r)rw(x.  v,  z  =  const)] 
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where  VV(  f\ .  f\  I  denotes  the  F  ourier  transform  of  h  i  v.  v.  r  - 
const)  and 

:  (I?) 

with  the  linear  frequencies  /,  and./,  corresponding  to  the  r 
and  v  directions.  F)y  treating  each  gA  as  a  planar  harmonic 
function  one  has  /.,igAt  -  agA  o:  ;md  /.,(gA)  -=  ( I 
r!)a’(  f>k  t  <i-' ,  Inserting  (7)  into  (b)  \  ields  (in  light  of  (12)) 


get  back  up  (or  back  down)  to  P.  One  gets  back  to  Ag /•  from 
Ag).  via  ( I )  w  hich  now  can  be  w  ritten  as 
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where  (/*'  denotes  the  two-dimensional  Fourier  transform  of 


Each  //  '  I  term  in  the  first  //-indexed  summation  is  obtained 
from  ( 14)  (recall,  by  definition,  g}t  =  Ag/>).  By  treating  each 
gA  as  harmonic,  it  follows  a  r(gA  )/<);'  =  F  '{G'A(-27r(/)r}, 
and  through  routine  algebra  one  can  easily  show  that  the 
second  //-indexed  summation  in  (21).  which  has  a  starting 
index  of  zero,  is  equivalent  to  the  following  sum  having  a 
starting  index  of  one: 


Fhe  well-known  analytically  defined  spectra  of  the  geo¬ 
metrical  kernels  of  (II)  are  given  by 


where  each 
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for  q  /  0.  Fort/  -  0  both  spectra  are  set  equal  to  zero  which 
corresponds  to  removing  the  mean  of  the  kernel. 

1. inear  convolution  in  the  space  domain  corresponds  to 
multiplication  in  the  frequency  domain  provided  the  neces¬ 
sary  steps  are  taken  to  convert  periodic  or  circular  convo¬ 
lutions.  resulting  from  the  use  of  the  direct  and  inverse  FFT. 
to  the  desired  linear  forms  (see  the  next  section  for  details). 
Hence  the  // th  linear  convolutions  of  ( 1 1 )  can  be  immediately 
expressed  (in  radians)  as 


Note  that  obviously,  for  «>  I,  g",  =  -g/i  and  (21)  reduces 
to  the  defined  A.g/>  =  gjl-. 

Similarly,  if  one  treats  each  £k.  £k.  and  17 A  as  planar 
harmonic  functions,  then  it  follows  that 

rT(c k)li>zr  =  [|/(27ry)]F  '{( l/</)(7*(  —  27rz/)r} 

=  ( -  I / y ) f '  '{Gk(-2irqY  -  '} 
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F  or  a  point  A  that  resides  on  both  the  Earth's  surface  and  the 
level  surface  (see  F  igure  I)  the  desired  Earth  surface  deflec¬ 
tions  are  simply  given  by 


f/r( T)k)Ut:r  =  (l ly)F  '(< —  (7,  lq)Gki -2nq)r\ 

=  t\!y)F~'{Gkt-2nq)r  '(2nift)\ 

Taylor  series  akin  to  (21 )  can  be  expanded  for  £P.  £P.  and  pr. 
Equivalent  to  these  expansions  are  the  following  final  deflec¬ 
tion  expressions  at  any  point  P  on  the  Earth's  surface: 


with  each  £'{  and  17"  computed  via  ( 17).  It  is  well-known  that 
the  analy  tical  Fourier  transform  of  |.v~  »  v’|  l_  is  l/r/forr/ 
'  0.  and  it  is  again  set  to  zero  for  </  -  0.  Hence  the  //th 
convolution  of  (10).  evaluated  (in  meters)  at  an  arbitrary  Q’ 
on  the  level  surface,  can  be  written  as 
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and  at  A  the  desired  Earth  surface  height  anomaly  is  given  by 


and  the  terms  in  the  first  //-indexed  sum  are  obtained  from 
(17).  For  the  height  anomaly  one  has 


i,.=  x  a  ‘  x  a 


II  n  I 


For  arbitrary  points  P  on  the  Earth's  surface  located  below 
or  above  the  level  surface  one  can  perform  a  Taylor  series  to 


'y  C; 

f  5:7 
3  gc 


■  or 

ziat  |  Special 


□  □ 


67X2 


Gi  fason:  Obtaining  Hakth  Surface  and  Spatiai  Dfi  i  potions 


f  ig.  2  The  1  arc  nun  by  I  arc  min  topography  (Oklahoma)  of 
central  1.75'  by  2.25'  area. 
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and  the  terms  in  the  first  /(-indexed  sum  are  obtained  from 
(19).  The  fact  that  the  efficient  and  simultaneous  computa¬ 
tion  of  the  Gk  spectra  only  has  to  be  executed  once  (the  Gk 
can  be  stored  on  a  direct  or  random  access  disc  file  for  future 
reference),  together  with  the  simultaneous  calculation  of 
gridded  £.  £.  and  17  values  via  (23M24).  constitute  the  appeal 
of  Sideris'  approach. 

3.  Potential  Pitfalls  of  the  Discrete  FFT  Pair 

Although  some  nice  treatments  of  geodetic  FFT  tech¬ 
niques  already  exist  [cf.  Schwarz  et  al .,  1990],  the  purpose  of 
this  section  is  to  furnish  a  short  (but  detailed)  explanation  on 
the  optimal  computer  execution  of  equations  (23M24)  using 
a  direct  or  random  access  file  of  precomputed  Gk  spectra 
values  obtained  recursively  from  F  (left-hand  side  of  (14)}. 
First  of  all.  most  "canned"  two-dimensional  FFT  routines 
such  as  International  Mathematical  and  Statistical  Library’s 
(IMSL)  FFT  3D  expect  the  Ag(0,  0)  origin  value  to  be  at  the 
southwest  corner  of  the  grid.  Thus,  to  make  inverse  Fourier 
transforms  such  as  (14),  (23),  and  (24)  amenable  to  such 
routines,  one  must  apply  the  appropriate  frequency  shifts  in 
any  analytically-defined  transfer  function  contributing  to  the 
final  spectra  which  will  be  transformed  back  to  the  space 
domain  by  the  inverse  two-dimensional  FFT  routine. 

There  are  three  major  areas  of  concern  in  any  use  of  the 
discrete  two-dimensional  FFT  pair,  namely,  aliasing,  spec¬ 
tral  leakage,  and  circular  (nonlinear)  convolution  effects. 
Aliasing  can  only  be  truly  eliminated  by  sampling  the  signal 
at  a  rate  at  least  twice  as  high  as  the  highest  frequency 
present  in  the  signal  (if  known).  Thus  one  can  only  hope  to 
minimize  the  aliasing  effects  by  using  as  dense  a  grid  of  data 
as  practically  possible.  A  popular  remedy  [cf.  Bergland 
1969]  for  spectral  leakage,  resulting  from  the  failure  to 
examine  the  signal  over  its  entire  assumed  period,  is  used 
here;  namely,  multiplying  the  given  set  of  gridded  free-air 
gravity  anomalies  by  the  well-known  two-dimensional  80% 
cosine  taper  window.  This  eliminates  discontinuities  along 
the  edges  of  the  grid  and  allows  the  data  to  take  on  some 


semblance  of  periodicity.  Such  a  window  must  be  applied  to 
the  actual  data  and  not  to  any  artificial  data  generated  by 
filling  out  a  record  with  zeros.  Such  artificial  data  will  be 
discussed  shortly  and  will  be  used  in  the  test  cases  of  section 
4. 

When  the  desired  linear  convolutions  of  (9) — ( 1 1)  are  com¬ 
puted  with  the  aid  of  the  discrete  periodic  FFT  pair,  such  as 
in  (14),  (23),  and  (24),  the  resulting  convolution  is  periodic  or 
circular  in  two  dimensions.  To  obtain  the  desired  linear 
convolution,  Oppenheim  and  Schafer  [1975]  show  the  di¬ 
mensions  of  the  two  sequences  to  be  convolved  must  be 
chosen  large  enough  to  avoid  unwanted  aliasing  inherent  to 
the  periodic  mode  of  the  FFT  pair.  Convolution  of  two  M  by 
N  sequences  produces  a  sequence  of  area  2M-I  by  2N-\. 
Hence  each  array  should  be  dimensioned  2 M  by  IN  (or 
larger).  Zero  values  are  assigned  to  each  location  in  the 
extended  areas  of  the  space  domain  gk  array.  If  one  evalu¬ 
ates  the  spectra  of  the  geometrical  kernel  or  response 
function  via  a  (southwest  origin)  two-dimensional  discrete 
FFT,  a  “four  corner"  space  domain  grid  assignment  is  made 
[see  Oppenheim  and  Schafer,  1975]  to  avoid  the  “wrap¬ 
around"  or  circular  convolution.  Therein,  the  corner  Nf. 2  by 
Ml 2  nonzero  sections  of  the  overall  2 N  by  2 M  grid  contain 
the  appropriately  shifted  space  domain  kernel  values  and  the 
rest  of  the  grid  is  zero.  One  then  takes  the  direct  FFT  of  the 
extended  2M  by  2N  (or  larger)  grids.  As  alluded  to  earlier, 
somewhat  better  predictions  occurred  by  using  analytically- 
defined  kernel  spectra.  Thus  the  analytically-defined  transfer 
functions  appearing  in  (14),  (23),  and  (24)  require  “extend¬ 
ed"  frequency  counter  shifts  (recall  the  transfer  function 
frequencies  must  be  shifted  to  adjust  for  the  southwest 
origin).  An  inverse  FFT  of  the  extended  product  spectra  is 
taken  and  the  desired  linear  convolution  is  located  in  the 
original  (nonextended)  grid  locations.  Experiments  have 
shown  that  neglecting  the  circular  convolution  effects  can 
increase  the  error  in  individual  deflection  predictions  by  as 
much  as  100%  and  increase  the  overall  RMS  of  the  predic¬ 
tion  errors  by  several  tenths  of  an  arc  second.  An  additional 
benefit  of  extending  the  original  space  domain  gk  grid  with 
zeros,  to  form  a  2 M  by  2 N  grid,  is  the  reduction  of  unwanted 
“picket  fence”  effects  in  the  Gk  spectra  (see  Bergland  [  1 969] 
for  more  details).  These  unwanted  effects,  namely,  the 
Fourier  coefficients  taking  on  the  behavior  of  a  set  of 
overlapping  band-pass  filters,  occur  because  the  gk  signal 
does  not  consist  of  a  desired  set  of  purely  discrete  orthogo¬ 
nal  frequencies. 

4.  Test  Results 

Sideris’  two-step  continuation  scheme  was  subjected  to  all 
of  the  procedural  steps  outlined  in  section  3  to  compute 
deflection  of  the  vertical  components  from  I,  3,  and  5  arc 
min  gridded  data  bases  of  mean  free-air  gravity  anomalies 
and  mean  terrain  elevations  in  two  areas  having  dramatically 
different  topographies.  The  I  arc  min  gridded  files  were 
furnished  by  the  Defense  Mapping  Agency  Aerospace  Cen¬ 
ter  (DMAAC),  St.  Louis,  Missouri.  The  coarser  3  and  5  arc 
min  files  used  were  simply  averaged  versions  of  the  I  arc  min 
data. 

The  north-south  Ay  planar  sampling  spacing,  in  meters, 
was  set  to  an  arc  length  of  I.  3,  and  5  arc  min  along  a 
meridian  of  radius  6371  km.  The  east- west  A.v  was  set  to  Ay 
•  cos  d»M.  where  <b M  is  the  middle-latitude  value  of  the 
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J-ijj .  3.  I  he  I  lire  min  by  I  me  min  topography  ( New  Mexico)  of 
central  3  by  V  area. 


surveyed  area.  Predicted  values  w ere  then  compared  to  a  set 
of  astrogeodetic  deflection  "truth"  values  located  in  the 
center  region  of  the  surveyed  area  (noncentral  predictions  of 
course  suffer  from  the  shortened  extent  of  the  data  and  from 
the  periodic  nature  of  the  FFT  pair  (with  respect  to  the 
hitter,  see  Warn;  1 1 9X8 1 ) . 

The  two  areas  selected  were  (I)  the  topographically  tran¬ 
quil  hut  gravimetrically  turbulent  (due  to  the  crustal  density 
structure)  Oklahoma/Texas  test  area  of  the  Gravity  Gradi- 
ometer  Survey  System  and  (2)  the  highly  rugged  White 
Sands,  New  Mexico,  area.  The  input  data  bases  for  Oklaho¬ 
ma Texas  covered  a  3.5  (north-south)  by  4.5°  (east-west) 
area  and  for  New  Mexico  a  6J  square  area.  Figure  2 
illustrates  the  central  1.75  by  2.25"  region  of  the  Oklahoma 
site  and  Figure  3  illustrates  the  central  3  by  3  area  of  the 
New  Mexico  site  (all  astronomically  determined  “truth" 
stations  were  located  in  these  central  regions).  Due  to  the 
rugged  nature  of  the  terrain,  the  accuracy  of  the  New 
Mexico  gravity  anomalies  is  far  from  uniform,  li  inaccessi¬ 
ble  areas,  a  least  squares  collocation  scheme  was  employed 
by  DMAAC  to  obtain  interpolated  anomaly  values. 

The  long-w  avelength  features  of  the  Farth  surface  gravity 
anomalies  corresponding  to  harmonic  degrees  n  ■  NT  were 
removed  by  the  spherical  harmonic  expansion 
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I II 


•  X  IT*,,,  cos  ill  A  t- .S„  sin  ///A  )P„  „,(cos  H)  (25) 

rtl  0 

where  H  denotes  colatitude  and  the  asterisk  implies  that  the 
potential  coefficients  were  made  residual  to  the  reference 
ellipsoid  (various  potential  coefficient  sets  and  values  for  NT 
were  selected;  the  results  will  be  discussed  below).  The 
geodetic  ( cf> .  A .  height)/,  coordinates  are  transformed  to  (r.  H, 
A)/,  geocentric  coordinates  before  (25)  is  executed.  The 
reduced  set  AgrcJllccd  =  A.g  -  Agrcmmcd  was  then  subjected 
to  the  two  step  continuation  scheme,  producing  a  set  of 
reduced  deflection  components.  The  long-wavelength  fea¬ 
tures  of  the  deflections  were  then  added  back  via 
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For  both  areas  the  mean  elevation  surface  was  chosen  as  the 
level  reference  surface  containing  the  points  P’  and  Q‘  on 
Figure  I  (experiments  revealed  the  /(-indexed  series  of  (5). 
(23).  and  (24)  converge  fastest  for  such  a  level  surface 
choice).  The  final  values  £t  and  t)a-  to  be  compared  to  the  Ath 
astrogeodetic  "truth"  values,  were  obtained  from  a  simple 
interpolation  of  the  four  predicted  (gridded)  deflections 
closest  to  the  Ath  control  point  given  by 


V  (I//),*) 


■  i 


where  D,t  is  the  planar  distance  from  the  control  point  A  to 
the  grid  point  i. 


I  Mil  I.  I  Maximum  e".  f'l  £/•  •  £/>).  and  ij" I  >;/”  •  rj/ll  Absolute  Values.  Given  by  (14)  and  (23).  in  Central  Region  of  Surveyed 

Oklahoma/Texas  Area 
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TABI.H  -  Overall  RMS  Values  of  24  (True  Predicted)  t  and  ?i  Differences  tor  the  Oklahoma/ 

Texas  Area 


n  Truncation 
Level  of  (23) 

n 

Grid  Resolution 

5  are  min  by  5  are 
min 

3  are  min  by  3  are 
min 

I  are 

min  by 
min 

1  are 

£ 

n 

£ 

V 

£ 

n 

0 

0.56 

0.42 

0.44 

0.42 

0.35 

0.32 

1 

0.55 

0.41 

0.43 

0.42 

0.34 

0.32 

■> 

0.55 

0.41 

0.43 

0.42 

0.34 

0.32 

3 

0.55 

0.41 

0.43 

0.42 

0.34 

0.32 

The  £  and  )j 

differences  are  in 

arc  seconds. 

Note  that  the  £,  "analytically  continued"  deflections 
represent  the  deviation  of  the  actual  plumb  line  from  the 
normal  plumb  line  at  the  ground  point  P.  while  the  astrogeo- 
detic  £  values  denote  its  deviation  from  the  perpendicular  to 
the  ellipsoidal  footpoint  of  P.  The  difference  is  accounted  for 
by  the  well-known  normal  reduction  for  the  curvature  of  the 
plumb  line  |cf.  Hciskanen  and  Moritz,  |%7|.  One  has 
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where  ./*  denotes  the  gravimetric  flattening,  R  the  radius  of 
the  spherical  Earth,  It  the  elevation  of  P,  and  6  its  geodetic 
latitude. 

For  the  Oklahoma/Texas  area,  use  of  the  Rupp  [  1 98 1 1 
potential  coefficient  set  through  NT  -  36  in  equations 
(25M27)  yielded  the  best  quality  predictions.  Although  such 
a  truncation  implies  a  5°  by  5°  local  data  set  is  available 
(recall  a  3.5°  by  4.5  set  was  used),  using  higher  degree 
coefficients  was  detrimental  to  the  final  predictions  ( Sidcris 
|I987]  cites  similar  experiences).  Current  high  degree  and 
order  geopotential  coefficient  sets  (such  as  Rapp's)  are 
"combined”  sets  in  the  sense  that  the  lower  degree  (n  s  36) 
coefficients  are  related  to  satellite  observations,  while  the 
higher  degree  coefficients  come  from  a  global  set  of  1°  or  30 
arc  min  surface  mean  gravity  anomalies.  The  latter  coeffi¬ 
cients  are  not  as  reliable  as  the  former  due  to  a  scarcity  of 
surface  gravity  information  over  large  portions  of  the  globe. 
Experiments  revealed  that  use  of  dense,  localized  gravity 


anomaly  data  sets  in  conjunction  with  the  satellite-related  n 
<  36  potential  coefficients  led  to  the  best  results. 

Table  I  lists,  for  the  5.  3.  and  I  arc  mir.  giids.  the 
maximum  ,g".  £"  (-  r  £?,)  and  17"  (-  17 *  rj/'i).  n  -  0. 
1.  2.  3.  absolute  values  in  the  central  region  of  the  Oklahoma 
area.  For  such  a  tranquil  topography,  the  higher-order  terms 
are  expectedly  impotent  (the  significant  g 1 .  £ 1 .  and  17 1  values 
(for  the  1  arc  min  grid)  were  mainly  located  near  the 
protruding  hills  appearing  in  the  middle  of  Figure  2.  The 
overall  RMS  (in  arc  seconds)  of  24  (true  -  predicted) 
deflection  differences  in  the  Oklahoma/Texas  area,  for  vari¬ 
ous  grid  resolutions  and  truncation  levels  of  the  /(-indexed 
sums  of  equation  (23).  appear  in  Table  2.  As  expected,  the 
higher-order  terms  contributed  little  to  the  final  interpolated 
deflection  computations.  The  average  standard  error  of  the 
astrogcodetic  "truth"  values  was  exceptionally  low,  =0.05 
arc  sec.  The  average  standard  errors  of  the  input  1  arc  min 
mean  anomalies  and  elevations  were  3.7  mGal  and  30  m. 
respectively. 

For  the  New  Mexico  site,  the  Rapp  [1981]  coefficient  set 
was  again  used  with  NT  =  36  to  account  for  the  long- 
wavelength  contributions.  Table  3  is  the  New  Mexico  coun¬ 
terpart  to  Table  i.  The  rugged  New  Mexico  topography  of 
Figure  3  is  reflected  in  the  large  magnitudes  of  the  Table  3 
higher-order  maximum  absolute  values.  Extrapolating  the  1 
arc  min  by  1  arc  min  section  of  Table  3  (it  was  truncated  at 
n  =  8  due  to  a  shortage  of  disk  space  at  the  Geophysics 
Laboratory  (GL))  infers  that  equations  (5)  and  (23)  should  be 


TABLE  3.  Maximum  g” .  f"(=fp+  <£).  and  ij"  (-Tjjl  +  rj/i)  Absolute  Values.  Given  by  (14)  and  (23),  in  Central  Region  of  Surveyed 

New  Mexico  Area 
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85.79 

2.19 

3.51 

■> 

3.31 
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0.16 

0.60 
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0.00 
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0.00 
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0.41 
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0.00 
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Fig.  4  The  e1  values  of  central  5°  x  V  New  Mexico  area, 
corresponding  to  I  arc  min  x  1  arc  min  grid. 


continued  through  n  =  10  to  obtain  submilligal  accuracies  for 
the  entire  level  surface  Ag'  set  and  0.1  arc  sec  accuracies  for 
the  Faith  surface  t  and  t)  sets.  Such  a  truncation  is  far  higher 
than  reported  by  Sideris  [  1 987 )  using  5  arc  min  by  5  arc  min 
gridded  data  bases  and  interpolated  I  arc  min  grids.  Wang 
1 1988|  also  reports  on  the  inability  of  interpolated  grids  to 
account  for  higher-frequency  terrain  effects.  Figures  4,  5, 
and  6  plot  the  g 1 .  g:.  and  g1  values  (vertical  scale  exagger¬ 
ated)  corresponding  to  the  I  arc  min  by  1  arc  min  grid,  for  the 
central  .V1  by  3  area  of  the  New  Mexico  site. 

Table  4  is  the  New  Mexico  counterpart  to  Table  2.  As  a 
preliminary  remark,  it  should  be  mentioned  that  the  378  pairs 
of  astrogeodetic  d  flections  used  are  not  part  of  an  older 
astro  set  used  by  the  well-known  C.  C.  Tscherning  White 
Sands  Study  Group  and  supplied  by  DMA's  Geodetic  Sur¬ 
vey  Squadron  in  Cheyenne.  Wyoming.  My  initial  experi¬ 
ments  using  the  Tscherning  astro  values  detected  a  bias  of 
approximately  1.2  arc  sec  in  the  east-west  tj  astro  compo¬ 
nents.  The  Survey  Squadron  subsequently  supplied  me  with 
a  newer  astro  set  having  accuracies  comparable  to  those  of 
the  Oklahoma  set.  The  values  appearing  in  Table  4  refer  to 
this  newer  bias-free  set.  The  tabular  results  agree  closely 
with  those  obtained  using  the  older  astro  set  after  the 
aforementioned  east- west  bias  was  removed. 

Although  the  modest  overall  improvement  in  the  predic¬ 
tions  from  the  use  of  the  I  arc  min  by  I  arc  min  g1  terms  is 


Fig  s  The  ir  values  of  central  3  >  y  New  Mexico  area, 

corresponding  to  I  arc  mm  «  I  art  min  grid. 


Fig.  6.  The  g1  values  of  central  3°  x  3°  New  Mexico  area, 
corresponding  to  I  arc  min  x  I  arc  min  grid. 


disappointing,  it  can  readily  be  explained.  Recall  gp  = 
~zpidgulr)z)p.  As  stated  earlier,  the  majority  of  the  total 
value  of  the  anomalous  vertical  gradient  is  due  to  the 
behavior  of  g  in  the  immediate  3  arc  min  area  surrounding  P. 
A  I  arc  min  grid  supplies  nine,  I  arc  min  by  I  arc  min  g° 
values  in  this  area.  For  near-mountain  astro  stations  having 
a  dense  sampling  of  actual  point  gravity  measurements  in 
this  critical  immediate  3  arc  min  area,  the  g1  term  consis¬ 
tently  enhanced  the  predictions  by  several  tenths  of  an  arc 
second.  By  contrast,  for  near-mountain  stations  lacking  such 
measurements,  the  g 1  term  tended  to  deteriorate  the  predic¬ 
tions  by  comparable  amounts.  The  overall  statistic  hides  this 
significant  g 1  activity  and  the  use  of  interpolated  anomalies 
in  areas  immediately  surrounding  the  computation  point  is 
detrimental  to  the  cause.  Harrison  and  Dickinson  [1989] 
suggest  constructing  a  30  arc  sec  grid  (in  this  case)  and 
placing  the  pertinent  I  arc  min  gravity  and  height  values  in 
each  of  the  four  30  arc  sec  cells  before  computing  the 
deflections.  Execution  of  this  suggestion  with  the  given  1  arc 
min  data  sets  did  rv.i  improve  the  predictions. 

In  light  of  the  above,  one  must  question  the  accuracy  of  g~ 
and  higher  terms  obtained  from  a  1  arc  min  data  grid.  The  g2 
term  can  justifiably  be  described  as  the  "gradient  of  the 
gradient"  of  anomalous  gravity.  Clearly,  the  majority  of  its 
value  is  accounted  for  by  the  behavior  of  g  in  an  immediate 
area  far  smaller  than  3  arc  min.  Higher-order  g*  terms  are 
even  more  short  wavelengthed.  One  must  question  how  an 
initial  set  of  I  arc  min  mean  anomalies  and  heights  can 
account  for  all  the  high-frequency  information  in  the  to- 
be-sol  ved-for  g*.  even  if  the  I  arc  min  mean  values  are 
straight  averages  of  a  dense  set  of  point  measurements  (i.e.. 
they  are  interpolation-free).  Moreover,  noise  (errors) 
present  in  the  I  arc  min  data  will  propagate  into  the  com¬ 
puted.  weak-signalled  gk  and.  in  all  probability,  render  them 
meaningless  due  to  unacceptable  signal/noise  ratios.  A  30  arc 
sec  grid  might  allow  for  meaningful  g~  values  provided 
interpolations  are  avoided  (of  course  one  runs  into  numerical 
stability,  computer  storage  and  CPU  concerns  with  such  a 
dense  grid).  By  the  same  argument,  if  one  is  working  with  3 
or  5  arc  min  data  grids,  one  is  hard  pressed  to  justify  even  the 
use  of  g1.  Notice  in  Table  4  that  the  overall  RMS  of  the 
east-west  rj  differences  actually  inched  upward  with  the  use 
of  the  3  and  5  arc  min  g1  values. 

As  a  final  remark,  one  can  easily  obtain  a  gridded  set  of 
deflections  at  a  constant  altitude  li  above  the  chosen  level 
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TABLE  4.  Overall  RMS  Values  of  378  (True  -  Predicted)  £  and  i;  Differences  for  the  New 

Mexico  Area 


n  Truncation 

Grid  Resolution 

5  arc  min  by  5  arc 
min 

3  arc  min  by  ?  arc 
min 

1  arc  min  by  1  arc 
min 

Level  Ol  1  * 

n 

i 

V 

e 

V 

€ 

V 

0 

0.84 

0.97 

0.72 

0.74 

0.61 

0.72 

1 

0.83 

0.98 

0.70 

0.76 

0.59 

0.66 

2 

0.83 

0.98 

0.70 

0.76 

0.59 

0.65 

3 

0.83 

0.98 

0.70 

0.76 

0.59 

0.65 

The  4  and  rj  differences  are  in  arc  seconds. 


surface  through  the  routine  spectral  execution  of  the  planar 
upward  continuation  integral  given  by 

{flat  height  h)\  =  F~  '{XAe  ~  2”hq} 

and 

{17 (at  height  h)}  =  F~  '{HAe  ~ 

where  XA  and  HA  denote  the  two-dimensional  Fourier 
transforms  of  the  level  surface  £A  and  tja  gridded  sets  given 
by  (18)  and  e~2rrhq  is  the  well-known  upward  continuation 
transfer  function. 

5.  Summary 

An  efficient,  spectrally  implemented,  density-free  tech¬ 
nique  allows  for  the  prediction  of  deflections  of  the  vertical 
and  height  anomalies  from  input  grids  of  free-air  gravity 
anomalies  and  elevations.  Long-wavelength  contributions 
are  accounted  for  by  a  low  degree  set  of  geopotential 
coefficients.  The  scheme  has  been  successfully  tested  in  two 
diverse  topographies  and  accounts  for  severe  terrain  effects 
as  well  as  turbulent  gravimetric  behavior  due  to  crustal 
density  structure.  One  minute  grids  of  noninterpoiated  free- 
air  gravity  anomalies  and  heights  are  required  for  a  universal 
set  of  accurate  first-order  vertical  gradients  of  anomalous 
gravity,  intermediate  quantities  in  the  prediction  process. 
Accurate  higher-order  gradients  will  require  even  finer  res¬ 
olution  grids  of  a  noninterpoiated  nature 
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